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S. V. FallertE 

Department of Chemistry, University of Cambridge, Cambridge, UK 

Y. M. Kim 

St. Catharine's College, University of Cambridge, Cambridge, UK 



oo 

o 

O 

(N 

C 

o 

S 1 

o 

a 
a 

C 

o 

> 

cn 

o 
l> 
o 



X 



,5Y. 



C. J. Neugebauer 
Department of Chemistry, University of Cambridge, Cambridge, UK 

S. N. Taraskin 

Catharine's College and Department of Chemistry, University of Cambridge, Cambridge, UK 

(Dated: June 30, 2008) 



The critical behavior of the contact process in disordered and periodic binary 2d-lattices is inves- 
tigated numerically by means of Monte Carlo simulations as well as via an analytical approximation 
and standard mean field theory. Phase-separation lines calculated numerically are found to agree 
well with analytical predictions around the homogeneous point. For the disordered case, values of 
static scaling exponents obtained via quasi-stationary simulations are found to change with disorder 
strength. In particular, the finite-size scaling exponent of the density of infected sites approaches a 
value consistent with the existence of an infinite-randomness fixed point as conjectured before for 
the 2d disordered CP. At the same time, both dynamical and static scaling exponents are found 
to coincide with the values established for the homogeneous case thus confirming that the contact 
process in a heterogeneous environment belongs to the directed percolation universality class. 



I. INTRODUCTION 

The contact process (CP) is a prototype model 
for the spatial spread of epidemics in biological systems. 
It describes epidemics in populations where each mem- 
ber can be in one of two states: infected (I) or sus- 
ceptible (S) (so-called SIS models). The CP exhibits a 
non-equilibrium phase transition between an active and 
a non-active regime of the disease, behaving at its criti- 
cal point according to the directed percolation (DP) uni- 
versality class. This has been established by a range of 



analytical and numerical techniques 0, H, 0, Hi such as 
renormalization group analysis [H, series expansions 
Q, Monte Carlo (MC) simulations [§, B and spectral 
analysis of the Liouville operator [Tol. These analy- 
ses have been undertaken for simple topologies, mostly 
for homogeneous hyper-cubic lattices. 

Recently, interest has turned towards the behavior of 
this process in disordered environments and revealed very 
peculiar features such as changing exponents and signif- 
icantly different dynamics like Griffiths phases and acti- 
vated scaling [l^, 0, 01 ■ In general, heterogeneous envi- 
ronments are typical in realistic systems, especially in the 
context of control of epidemics [la, Us LLli ■ Therefore, 
it is instructive to investigate the critical behavior of the 
CP under these conditions and in particular to establish 
the phase diagrams for such systems. In the past, the dis- 
ordered CP (DCP) has been investigated in both one and 
two dimensions in a range of settings and revealed a con- 
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tinuous change in static critical exponents starting from 
the clean DP values [H, [H, H3] . I n the Id case, a strong- 
disorder renormalization group study allowed deep in- 
sight into the disordered process and revealed a domi- 
nating infinite-randomness fixed point (IRFP) for suffi- 
ciently strong disorder Q . In the weak- to intermediate- 
disorder regime however, MC simulations and density- 
matrix renormalization g roup (DMRG) techniques [l9| as 
well as series expansions [2fJ found continuously varying 
disorder-dependent critical exponents which were found 
to approach those characteristic of an IRFP with increas- 
ing strength of disorder. For the 2d CP with site dilution, 
MC simulations showed a similar behavior, a continuous 
change in exponents with increasing disorder [12l . I IL 31 ] and, 
in retrospect [19| , giving evidence for the existence of an 
IRFP also in this case. 

In this paper, we consider the phase diagram of the 
CP in a binary 2d lattice of sites with different recov- 
ery rates ca and cb drawn from a bimodal probability 
distribution. Extensive Monte Carlo (MC) simulations 
following [l2| are employed in order to locate the line of 
critical points in the space of recovery rates (ca, £s)- As 
such simulations of disordered systems are numerically 
intensive due to very long relaxation times, analytical 
approximations are vital to constrain the region of phase 
space which contains the line of critical points. Here, we 
analyze the results obtained from the mean field approx- 
imation and further propose an approximate expression 
for critical points guided by the structure of the Liouville 
operator which governs the time evolution of the CP. 

Also, the quasi-stationary (QS) simulation method [2l| 
is employed to investigate the static scaling behavior of 
the CP in this disordered system and to calculate the cor- 
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responding critical exponents. In particular, we study 
whether the process in this setting exhibits disorder- 
dependent changing exponents which cross over to values 
characteristic of an IRFP for sufficiently strong disorder 
as observed previously [H, [l9| . 

Following on, we investigate the behavior of the CP in 
a range of heterogeneous periodic lattices with different 
unit cells via MC simulations and test the validity of our 
analytical expression as well as standard mean field the- 
ory. The two analytical approaches, mean field and our 
alternative approximation, enable us to largely constrain 
the location of the critical points in both the disordered 
and the heterogeneous periodic lattices. Critical expo- 
nents are found to change continuously in the former case 
with increasing disorder and appear to approach the pre- 
dicted values characteristic of an IRFP while they remain 
constant at their DP values in the latter. 

The CP, its critical behavior and some of the theoret- 
ical foundations employed for its description are intro- 
duced in Sec. [TTJ Our analysis of the disordered system 
is presented in Sec. IIIII In Sec. [IV] we investigate the CP 
in a range of heterogeneous periodic lattices in a similar 
fashion. Lastly, our findings are discussed in Sec. [V] and 
we summarize in Sec. IVI1 

II. BACKGROUND 

In this section, we define the CP and give an overview 
of its critical behavior and the master-equation descrip- 
tion by means of the Liouville operator. The CP is a 
non-equilibrium stochastic process in which an infection 
spreads via nearest-neighbor contact from site i to j at a 
transmission rate Wi—,j . Recovery of site i is spontaneous 
and happens at a recovery rate In the thermodynamic 
limit, the ratio of these two rates is the control parameter 
of a second-order phase transition between a non-active 
phase where no infected sites remain as t — > oo and an 
active phase where the density of infected sites (order 
parameter) is non-zero as t — ► oo 0, 0, 0| • 

For the CP in a system of size N with sites i = 1 . . . N 
we denote the two possible states of site i as Si = 1 
(infected) or Sj = (susceptible). A microstate of the 
system, i.e. a snapshot of the infection states of all sites, 
can be defined as a vector S = (si, . . . ,sn) T and the 
probability of finding the system in a specific microstate 
at time t is denoted by P(S, t). Assuming the transition 
rates between microstates S and S' to be rs— »s' , the time 
evolution of this probability follows the master equation 
which expresses the conservation of probability flow, 

d t P(S, t)=J2 (rs^sP(S', t) - r s ^ s ,P(S, t)) , (1) 

S' 

where the transition rates rs^s' follow from the rules of 
the CP. The master equation can be recast in compact 
form by introduction of the Liouville operator C which 
acts on the probability state vector \P(t)), 

d t \P{t)) = t\P{t)) , (2) 



the components of which are the probabilities of find- 
ing a system of N sites in different states \a) at time t, 

= E CT M-PM}^)- Here > iW)} is the orthonormal 
basis diagonal in the occupation number representation 
3, 10]. The precise form of L is most readily expressed 
in terms of hard-core bosonic creation and annihilation 
operators acting on site i, aj and ai, respectively, 

L = Y1 \ ei i 1 ~ a t) a i + i 1 ~ a i) a i w ^ a ) a i '( 3 ) 

i \ 3 eNN(€) J 

where the first part destroys particles while the second 
part creates offspring The Liouville operator is non- 
Hermitian with matrix elements = (cr' |jC|ct) , which 
coincide with the transition rates from state a to state 
a' ^ a and t aa = - J2a'^a Ar'<r- 

In principle, Eq. @ can be solved by performing direct 
diagonalization of the 2^ x 2 N real sparse (for lattice 
topologies) non-symmetric Liouville matrix. Its formal 
solution can then be expressed as 

\P(t)) = 5^(^(0)) \ei) , (4) 

i 

where are the eigenvalues of the Liouville matrix with 
a complete set of eigenstates |e») . The trivial solution 
|eo) of the eigenproblem for the Liouville operator with 
A = corresponds to the absorbing state of the system. 
All other eigenvectors |e$) in finite systems have eigenval- 
ues with negative real parts and thus decay exponentially 
with time. In the thermodynamic limit (N — * oo), there 
is one eigenstate |ei) with corresponding eigenvalue, Ai, 
which is zero in the active and non-zero in the non-active 
phase. In a finite system, the value of Ai in the active 
(non-active) regime approaches a zero (non-zero) value 
with increasing N, thereby signaling the phase transi- 
tion. The exact location of the transition can be extrapo- 
lated using finite-size data for moderate system sizes (e.g. 
N < 16) from direct diagonalization or density-matrix 
renormalization group calculations 0, [T3| • 

III. DISORDERED SYSTEM 

In what follows, we investigate the behavior of the CP 
on a lattice of two types of site, A and B, characterized 
by different recovery rates ca and es, respectively. The 
recovery rate at site i, ej, is drawn from the bimodal 
distribution 

p(ei) = x 8(e l - e A ) + (1 - x) 5 fa - e B ) , (5) 

where x controls the relative concentration of A and B 
sites. The transmission rate, for simplicity, is the same 
for all possible links between nodes, Wi—,j — Wj^i = w. 
As a further simplification, the timescale is set up by 
choosing w = 1 jZ with Z being the number of nearest- 
neighbor links per node (Z = 4 for the topologies consid- 
ered here). 
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MC Simulation 




case of the disordered system considered above, the gov- 
erning equations for the mean concentrations of infected 
sites of type A and B, n A and tib respectively, are given 
by 

= -e A n A + y(l - x)n B {l ~ n A ) + ^xn A {l - n A ) 

where = wZ = 1 is the mean field critical value 
for the recovery rate in the homogeneous 2d square lat- 
tice. As usual for the mean field approximation in 
low-dimensional systems, e* significantly overestimates 
the true critical value for the homogeneous case, e A = 
eg. The locus of critical points in the parameter space 
(e^jCs) separating non-active and active phases can be 
easily found from the solution of Eqs. ([6]) in the steady- 
state regime giving 
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(7) 



FIG. 1: (Color online) The phase diagram for the CP on a 
disordered lattice with sites of recovery rates ea or eb drawn 
from the distribution Eq. (JSJ for the case x = 0.5 obtained 
from MC simulation (dots), mean field (upper dashed line) 
and the analytical expression Eq. (JS| (solid line) . Inset shows 
the deviation A(d) as defined in the text for both mean field 
(blue □) and Eq. © (red o). 

MC simulations are used to locate the critical point by 
starting from a single infection seed and averaging over 
10 6 realizations of the process each with a fresh realiza- 
tion of the disorder up to a maximum of 10 6 time steps. 
We consider the case of a; = 0.5 and a range of values 
for e A , aiming to find the corresponding critical cb for 
each. As previously observed [12j, very slow dynamics 
are encountered, an effect that increases with disorder 
strength, i.e. the distance between e A and the homoge- 
neous critical rate e c = 0.60653(3) ,3]- Following [T^], the 
criterion of asymptotic monotonic growth or decay in or- 
der to assess whether the process is super- or sub-critical 
for a particular choice of rates is employed. For clarity of 
presentation, rescaled recovery rates = e,/ e c are intro- 
duced in order for the homogeneous critical point to be 
conveniently located at e c = e A = e B = 1- The result- 
ing phase diagram, symmetric in e A and es, is shown in 

Fig.ru 



with Ei = £i/e* where (...) denotes an average over dis- 
order realizations. 

The resulting phase separation line is shown in Fig. [1] 
along with the MC data presented above. Note that due 
to rescaling, mean field and numerical results coincide 
by construction at the homogeneous critical point. In or- 
der to allow a quantitative comparison between numer- 
ical results and approximation, a measure of difference 
between prediction and the true value obtained by MC 
simulation is needed. As such a measure, we consider the 
shortest distance A(d) between the prediction curve and 
an MC data point a (shortest path) distance d away from 
the homogeneous point. This error quantity is suitable 
for quantitative analysis as it is a measure for the width 
of the region of uncertainty between the analytical pre- 
diction and the true critical line and will be symmetric 
about the homogeneous point for symmetric phase dia- 
grams. The inset of Fig. [T] shows A(d) for the mean field 
approximation (blue squares). 

As can be seen from the figure, the mean field ap- 
proximation provides an upper bound to a region that 
contains the phase separation line. While the deviation 
A(g?) is small in the vicinity of the homogeneous critical 
point (A < 0.01), it grows considerably as the degree of 
heterogeneity increases (A w 0.1). 



B. Mean Field Approximation 

As outlined in the introduction, we are interested in 
analytically approximating the region where the phase 
separation line is located. To this end,_we first present an 
approach based on mean field theory [3J] . In this approxi- 
mation, fluctuations and correlations are ignored render- 
ing the master equation analytically tractable. For the 



C. Alternative Analytical Approximation 

Given that the mean field approximation appears to 
provide an upper bound to the region which contains the 
phase separation line, we are interested in obtaining an 
alternative analytical approximation that may provide a 
lower bound. In the following we will first present an ap- 
proximate expression for the location of critical points in 
a heterogeneous system and then compare its predictions 
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to the MC data of section UlI Al Following on, we give a 
motivation for this approximation along with numerical 
support. 

Statement and Comparison to Data 

Consider a finite system of N sites with arbitrary re- 
covery rates et (i = 1, . . . , N). We will argue below that 
for such a system in the vicinity of the homogeneous crit- 
ical point (all Ci = e c ) the expression 

JV 

rb= e *> (8) 

3 

approximately predicts the location of critical points. For 
the disordered system presented earlier, this expression 
simplifies to = ea £b- Fig. Q] shows both the MC 
data presented earlier and the approximate line of critical 
points thus obtained. As can be seen from the figure, the 
alternative analytical approximation is found to provide a 
reliable lower bound to the region which contains the line 
of critical points. Hence, in combination with the mean 
field approximation discussed above, one can constrain 
this region. Considering the error A(ii), it is found to 
show similar behavior to the one previously observed for 
the mean field data albeit an order of magnitude larger. 



Motivation and Numerical Support 

In the following we motivate Eq. ([5]) by considering the 
structure of the Liouville operator as defined in Sec. [Til 
For the finite system introduced above, the eigenvalues of 
the Liouville operator, A, are given by the characteristic 
equation, 

\- 1 \\I-£\=Q N _ x ({e l },\) = J2 A n ({ei})X n =0 , 

(9) 

where Q7v max (A) is a polynomial in A of order A^ max 
(iVmax = 2-^ — 1) and division by A eliminates the triv- 
ial zero root for the absorbing state. It is our aim to 
solve this equation approximately in the vicinity of the 
homogeneous critical point where = e c . 

To this end, we first consider the coefficients A n and 
look for features in their structure which may help in 
rendering the equation tractable. Generally, the A n can 
be expressed as 

N 

A n ({ei})= J2 'C IK • (10) 

mi r ..,mjY-0 j 

where the upper limits in the sum depend on n but their 
precise values are not significant for the analysis below. 
We now assume that in the construction of the A n ({ei}) 



from the determinant of XI — £, the dominant contribu- 
tion stems from terms with products of the same (or at 
least similar) powers of recovery rates at different sites. 
If this is true, the previous equation can be approximated 

as 



N \ 

IM - (ii) 

3 J 



where m* < 2 N /N. While this assumption may at first 
appear artificial, a justification can be found in the struc- 
ture of the Liouville operator. The determinant of the 
Liouville matrix contains the sum of terms which are the 
products of the recovery rates (the transmission rates 
are chosen to be constant). A typical (representative) 
term contains a product of many recovery rates, each 
one picked from a different column. Assuming periodic 
boundary conditions, all sites in the system should en- 
ter the Liouville matrix in the same fashion. Therefore, 
in a typical term one would not expect to find an over- 
representation of a specific site leading to the statement 
that the (combinatorially) dominant terms correspond to 
products of recovery rates raised to powers that are close 
in value. 

This argument only holds if one can be sure that the 
combinatorial weight of terms with homogeneous pow- 
ers is not offset by the actual values of the recovery rates 
&j. Otherwise, one could imagine the dominance of terms 
with very different powers of ej caused by the raising of 
values > 1 to a high power. However, recall that for 
the clean CP at the homogeneous critical point, the true 
critical recovery rate is < 1. Thus, close to this point, 
the critical recovery rates will always be close in value 
and < I which means that the above argument about 
homogeneous powers is expected to hold in this regime. 
Further support will be given below in the form of nu- 
merical evidence using a specific system further down. 

As explained in Sec[TTJ at criticality the highest non- 
trivial eigenvalue Xi(N) will be finite and tends to zero 
with increasing N. For the case of homogeneous recov- 
ery rates at criticality, e, = e c and Qiv maJt ({ e c}> Ai) = 0. 
Finally, by combination of this property and Eq. (1II[) . 
we indeed find Ylj e j — e c f° r t ne homogeneous critical 
point which is precisely the statement presented above in 
Eq. ©. 

The last step of our approximation then is to employ 
the same relation away from the homogeneous point and 
to use it to predict the locus of critical points. 

More formally, the above condition for critical points 
can be derived from Eqs. © and (TIT]) via a Taylor series 
expansion of Q./v max ({ei}, Ai) around the homogeneous 



5 



critical point in ln(ej/e c ), 

JV max 

Q;w(tehAi)= E Mia})*" 

- E E^ n)eTOEf H A ? 



n— \m— 1 



oo u 
171 



N 



E E^™"EiT^(nrl A "-^ 2 ) 



n— m— 1 



fc=l 



Here, we used the relation Qjv max ({ec}, Ai) = 
^ n ™Q X A„(e^)A" = which leads to no constant term 
in the expansion and allows factorization of the above 
expression, i.e. 

Qiwfeh Ar) = S In (j[ = , (13) 



where 



i v max iaj ^ / \ 

s = E E E lnfe " ( II ~ ) A " ■ ^ 



n— m— 1 



fe=l 



Eq. (O is obeyed if 



= , 



(15) 



because 5 ^ for arbitrary choice of which coincides 
with the condition given by Eq. ((5]). Alternatively, our 
approximate expression can be recast as an expectation 
value of logarithms [2(| , 



E 



ln-1 







(16) 



Note that this procedure amounts to simply geometri- 
cally averaging the recovery rates and inserting them into 
the clean theory. Interestingly, the logarithm of rates well 
known from renormalization group analyses of the DCP 
and the random transverse-field Ising model Q arises 
naturally in our scheme. 

In order to support the assumption about a dominant 
contribution to Eq. (| 10|) from products of homogeneous 
powers of recovery rates, numerical evidence for a simple 
system is given below. Let us consider a Id binary chain 
of sites A and B characterized by recovery rates 6a and 
es, respectively, and spatially arranged as . . . ABAB . . . 
with periodic boundary conditions. As a particular ex- 
ample, we analyze the coefficient A (eA,^B) defined by 
Eq. (JTDJ) , which reads (where a = a^) 

— a m A m B € B 

mA-.rriB 
n I * 



^ B m (e A ,e B ) 



(17) 



with 



B m (e A ,e B ) = (e A e B ) m x 



m* — m 



. (18) 



which can be symbolically evaluated for relatively small 
systems (N < 6). Initially, ea and e B will both be set 
equal to e c consistent with our assumption that we in- 
vestigate the vicinity of the homogeneous critical point. 
This enables us to investigate the relative magnitude of 
terms corresponding to different arrangements of powers. 
The terms B m effectively correspond to the contributions 
which contain either the homogeneous power m or one 
recovery rate to the power m along with the other recov- 
ery rate to a power greater than m. The magnitudes of 
the B m as functions of m for the binary system of size 
N — 4 and N — 6 are shown in Fig. [2] (top panel) . For 
both cases we observe sharp peaks centered at m max = 3 
(N = 4) and m max = 12 (N = 6) indicating a dominant 
contribution from a narrow range of powers. 



B 0.5 




k/m 



FIG. 2: (Color online) Top panel: The terms B, n as defined 
in Eq. (TT7)| for a linear AB chain of size N = 4 (left black 
peak, O) and N = 6 (right red peak, □) normalized by their 
maxima. Bottom panel: The correction terms Ck as defined 
in Eq. (|19[) as a function of the relative difference between 
powers k/m m ax where m ma x is the location of the respective 
maximum in the upper panel. Symbols as before, all val- 
ues have been normalized to the corresponding homogeneous 
contributions Co- 

The contribution to Aq from purely homogeneous pow- 
ers can be written as Co = J2m=o a m,m (e^es)™ while 
corrections to this can be expressed as 



Cfc — E ( a m+k,m £ 



m^m+k\ 
A e B ) 



(19) 
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for k > 1. The values of Ck represent contributions from 
powers differing by k from each other thus allowing a sys- 
tematic investigation of the validity of our assumption. 
We are interested in the magnitude of these corrections as 
a function of the relative difference normalized by m max 
between powers in order to allow a comparison between 
different system sizes. While the homogeneous contri- 
butions Co are found not to be the most dominant, the 
corrections are peaked at C\ for both systems considered 
and decay quickly with k. In particular, this decay hap- 
pens increasingly rapidly with larger N as a function of 
relative difference between powers, A;/m max (cf. the red 
curve marked by squares (□) for TV = 6 and the black 
one marked by diamonds (O) for N = 4 in Fig. [2]) (bot- 
tom panel) . A deviation of the values of ca and e b from 
their value of e c is found to reduce the dominance of the 
peaks presented above but does not immediately inval- 
idate the assumption. However, when moving far away 
from the homogeneous critical point, the peaks flatten 
out indicating a breakdown of our approximation. In 
summary, all of the above findings can be considered to 
support the assumption about a dominant contribution 
of homogeneous powers of recovery rates in Eq. (jTTJ) . We 
have undertaken a similar analysis for the coefficient Ai 
and expect the same behavior for the remaining A n . An 
analysis of A n (for n > 2) in a similar manner quickly be- 
comes prohibitive due to the computational complexity 
of the resulting expressions. However, as Ai approaches 
zero with increasing N, these higher terms are expected 
to become increasingly irrelevant. 

The question of whether one always expects to obtain a 
lower bound is addressed in the discussion (Sec. [Vj after 
more example cases have been compared to simulation 
data. 



D. Critical Exponents from Quasi-Stationary 
Simulations 

Investigations of the 2d DCP have been carried out in 
the past and have investigated both d yna mic [l2j and 
static scaling properties of the process [1J|. In general, 
the study of critical properties of the disordered process is 
complicated due to long relaxation times and ambiguity 
regarding the nature of scaling. In the following, we will 
investigate the static scaling of the disordered process by 
employing QS simulations [2l| and compare our results to 
both previous studies as well as theoretical predictions. 

In the clean CP, the order parameter, lim t _ (00 p is ex- 
pected to obey the scaling form Q 



(20) 



where x = f3/v±, L is the linear size of the system, and 
(3 and z/j_ are critical exponents. Further, G is a scaling 
function which asymptotically behaves as G{y) — > y@ as 
y — > oo and G(y) — > const, for y — > 0. An analogous 
finite-size scaling form is expected to be obeyed by the 



order parameter fluctuations, x — ^ d (^p 2 — p 2 ^j , with the 

exponent x replaced by x' — —*y/is±. 

In order to apply the above scaling relations, one com- 
monly considers QS values of observables as no true 
stationary state can exist in a finite system. The CP, 
when started from a fully infected system, initially re- 
laxes while spatial correlations grow towards the system 
size and temporal correlations decay. Once the spatial 
correlation length becomes comparable to the size of the 
system, the process enters a QS regime characterized by a 
time-independent non-zero transition rate to the absorb- 
ing state. In this regime, the QS density p, i.e. the den- 
sity p conditioned on survival, attains a constant value. 
In the past, analysis of this metastable state in com- 
puter simulations has proved to be notoriously difficult. 
Usually, the time-dependent density of infected sites con- 
ditioned on survival, p, which becomes stationary in the 
QS regime, is investigated Problematically though, it 
is neither clear at what time this density has converged 
to its QS value nor when the QS state starts to decay 
due to finite-size effects [22]. Therefore, a range of alter- 
native approaches have been proposed which enable an 
observation of this metastable regime (see Ref. [2l[ and 
references therein). Here, we employ the QS simulation 
method [2l| which allows a direct sampling of the QS 
state by eliminating the absorbing state and redistribut- 
ing its probability mass over the active states. 

Following Ref. [2l|, one starts from the master equa- 
tion Eq. (fTJ) . For the CP, this equation does not admit a 
non-trivial stationary solution for a finite system due to 
the existence of the absorbing state which can be en- 
tered but not be left. The QS solution mentioned above 
can be defined as 



P(S) = hm 

V ' t^oo P s (t) 



(21) 



where P s (t) denotes the survival probability of the pro- 
cess at time t. Now, consider a modification of the gov- 
erning equation, 

d t Q(S,t) = J2[rs'^sQ(S',t)-r s ^ s ,Q(S,t) 

S' 

+rs^ Q(S',t)Q(S,t)} . (22) 

where Q(S, t) denotes the probability of a new process 
governed by this equation being in state S at time t. 
The stationary solution of Eq. ([22]) . Q(S), coincides with 
the QS probability of the original process as can be seen 
by substituting Q(S,t) = P s (t)P(S) and noticing that 
in the QS regime dP s /dt = —P s ^ s rs^o-P(S). In that 
case, the right-hand side of Eq. (|22|) is equal to zero if 
Q(S) = P(S) as required. The last term in Eq. ([22]) 
can be viewed as a redistribution of probability from the 
absorbing state to the active states according to their 
probability [2lj . Thus, if one could sample from a pro- 
cess governed by Eq. (|22p . it would converge to a true 
stationary state governed by the QS probability distri- 
bution of the original process. Such a process is given 
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by the original CP where all transitions to the absorbing 
state are instead redirected to an active state randomly 
chosen according to its probability. As in practice this 
probability is not known a priori, an estimate is gener- 
ated by sampling from the history of the process. Gener- 
ally, the method has proved to be efficient with fast and 
reliable convergence after optimization of history sam- 
pling parameters [2ll [23J ■ The approach is particularly 
suited to a study of the DCP for which, in dynamic single- 
seed MC simulations employed for the DCP in the past 
[l2l [l4j , the question of whether the asymptotic limit of 
the process had been reached was frequently contested. 
In contrast, QS simulations offer a clear means of ensur- 
ing this: a true stationary average whose convergence can 
be monitored. 



TABLE I: The critical rates es for a given €a, critical ex- 
ponents x and x' for the disordered systems discussed in the 
text. 



points reads 





es (critical) 


X 


x' 


0.60653 


0.60653(3) 


0.795(4) 


0.42(3) 


0.595 


0.6188(3) 


0.796(5) 


0.41(5) 


0.5 


0.7676(4) 


0.83(1) 


0.39(3) 


0.4 


1.1815(5) 


0.92(4) 




0.35 


1.7775(5) 


0.93(4) 




0.3 


3.89(1) 


0.99(5) 





Here, we have investigated the 2d DCP with bimodal 
disorder in its recovery rates drawn from the distribu- 
tion Eq. ([5]) by means of QS simulations for up to 10 8 
time steps and systems of sizes L = 8, . . . , 128 sites av- 
eraging over no less than 10 3 disorder realizations. At 
the critical point, fits to the above finite-size scaling re- 
lations yielded estimates for the exponents x and x' . 
For the homogeneous case, the well-established values for 
the exponents of the DP universality class are recovered 
= 0.795(7), 7/V1 = 0.41(2) [H). As the degree 
of disorder, i.e. the difference between recovery rates ca 
and es, is increased, the measured exponents are found 
to change with disorder strength where x increases while 
x' decreases (cf. Tab. IJ). For strong heterogeneity, no 
credible fluctuation exponent could be extracted from the 
data due to strong sample-to-sample fluctuations. This 
is unfortunate as it prevents us from testing the validity 



of the hyperscaling relation 



2£ 



A similar re- 



lation for dynamical exponents had previously [12J been 
found to break down for the DCP. 



IV. 



HETEROGENEOUS PERIODIC LATTICES 



In order to investigate the range of validity of Eq. ((8]) , 
we now turn to the behavior of the CP on lattices with 
periodic arrangements of sites of type A and B discussed 
above. 

For such systems, the equation for the locus of critical 



(23) 



where ca and cb denote the concentration of species A 
and B, respectively. Three lattice systems have been an- 
alyzed with ca = cb = 1/2: (i) a standard chessboard 
lattice [Fig. 3(a)] , (ii) a big c hessbo ard lattice [Fig. 3(b)] 
and (iii) a lattice of rows [Fig. |3(c)| . Extensive MC simu- 
lations (3 x 10 6 runs up to t — 3000 maximum time steps) 
starting from a single infection seed were performed for 
these heterogeneous lattices. Unlike in the disordered 
case, for heterogeneous systems, asymptotic scaling re- 
lations that are well-known from the homogeneous case 
are found to hold. At criticality, the average number of 
infected sites (N(t)), the mean squared radius (R 2 ) of 
spread of the CP (where angular brackets denote aver- 
aging over all realizations and over active realizations at 
time t, respectively) and the survival probability P(t) 
follow asymptotic scaling laws [H, 



(N) ~ f> , (R 2 ) ~ t 



2/z 



(24) 



where 77, S and z are the dynamical critical exponents 
characteristic of the universality class. These scaling re- 
lationships provide a method for finding the critical value 
of the control parameter by fitting the observables to the 
above scaling forms following Ref. @. Furthermore, the 
dynamical critical exponents can be determined from the 
fit. 

As expected, the numerical data agree very well with 
the analytical predictions given by Eq. ([2"3"| [cf. the circles 
with the solid line for sa — £b — 1 in Figs. 3(a)|3(c) in 
the neighborhood of the homogeneous critical point, and 
start to deviate from the predicted phase-separation line 
for Sa,b 1 consistent with the validity of our approx- 
imation. The quality of the analytical approximation is 
high for the standard chessboard case (A < 0.03 for a 
very large range of rates) but becomes worse for the big 
chessboard (A up to 0.15 when moving away far from the 
homogeneous point) and especially for rows in the range 
of large values of ea.b S> 1. 

Furthermore, we have studied two lattices with differ- 
ent concentrations of nodes A and B, i.e. ca/cb = 2/1 
- lattice (iv) [see Fig. |3(d)| , and ca/cb = 3/1 - lattice 
(v) [see Fig. 3(e)] . In these cases, the phase-separation 
lines are not symmetric about the bisector in the sa — £b 
plane. As can be seen from Figs. 3(d)|3(c) the results of 
MC simulations of the CP on these lattices are again in 
good agreement with the analytical expression given by 
Eq. (|23p. especially near the homogeneous critical point 
(cf. the circles with the solid line for ea — £b — 1 in 
Figs. 3(d)||3(c)" ). In case of lattice (iv), the error as shown 
in the inset indicates a similar order of magnitude de- 
gree of accuracy as in the simple chessboard case before 
(A < 0.05 for a large range of rates) while for lattice 
(v) the approximation is found to deteriorate [with ap- 
proximately twice the value of A as compared to lattice 
(iv)]- 
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(d) Lattice (iv) (e) Lattice (v) (f) 



FIG. 3: (Color online) 

(a)-(e): The phase diagram for the CP on the lattices (i)-(v) as defined in the text. The circles represent the MC data, the solid 
line is given by Eq. (|23p for the specific lattice, and the dashed line corresponds to the mean field result taken from Tab. [TTJ 
Inset shows the deviation A(d) for both approximations (mean fields □, Eq. (I23[l = o). 

(f): The deviation A as defined in the text at criticality for e A = 0.9 for the CP on a lattice such as case (ii) but with variable 
linear size L of the square contiguous regions of A or B sites. 



TABLE II: The expressions for the phase-separation lines for 
the CP on different lattices (first column) obtained according 
to Eq. I|23[l (second column) and within the standard mean 
field approach (third column). 

mean field 



lattice type Eq. ([23 



(i) 
(ii) 
(iii) 
(iv) 

(v) 



£ B = 1/SA 

EB = 1/eA 

£ B = 1/£A 

£ B = 1/Ea 

e B = l/e A 



E B = 1/ea 

EB = £a/(2ea - 1) 
EB = £a/(2£a - 1) 
EB = 1/(2£A " 1) 
EB = E A /{2e 2 A - 1) 



It is instructive to compare the expression for the 
phase-separation line given by Eq. <P5|) with the re- 
sults obtained from the master equation within the stan- 
dard mean field approximation. Expressions similar to 
Eqs. © can be found for all the different lattices and 
solved for the critical rates in the steady-state regime. 
The resulting expressions for all the lattices are summa- 
rized in Table HIl 

As follows from Table UH the mean field result coincides 
with the expression for the phase-separation line given by 
Eq. P5|) for the standard chessboard configuration [lat- 
tice (i)] and gives a different prediction for all other cases 
studied. The rescaled mean field results agree very well 



with MC data around the homogeneous point but display 
deviations for ea.b ^ 1. Looking at the corresponding 
errors A(d), they are found to be of the same order as 
found for the previous analytical approximation. 

The fact that for the simple chessboard lattice our ear- 
lier prediction and the mean field result coincide reveals 
this case to be special in that the rescaled mean field does 
not over- but underestimate the true critical values. In all 
other studied lattices, the rescaled mean field results for 
the phase-separation lines lie above the numerical data 
[cf. the dashed lines with the circles in Figs. 3(b)||3(ej] and 
thus lead to an overestimate of the value if Eb for a given 
ea- This means that for these cases, mean field estimates 
of critical values can serve as an upper bound on the crit- 
ical recovery rate. In contrast, the phase-separation lines 
predicted by Eq. (J22J) provide a consistent underestimate 
of the true critical line for all studied lattices and there- 
fore a lower bound (cf. the solid lines in relation to the 



circles in Figs. 



and see the arguments given in 
Sec. [V} for the critical thresholds. 

In order to more systematically investigate how our 
alternative analytical approximation deteriorates as the 
spatial arrangement of sites becomes "less mixed", we 
consider a lattice like the big chessboard, lattice (ii), but 
vary the linear size L of the square contiguous regions 
of A or B sites. The resulting deviation A(d) as defined 
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above at critical cb for a choice of ea — 0.9, i.e. appre- 
ciably far away from the homogeneous critical point, as a 
function of L is shown in Fig. |3(f)| One can see that for 
L > 4 the accuracy quickly becomes worse than for any 
rate and lattice previously considered indicating a rapid 
breakdown of the approximation. 

Finally, the universality of the critical behavior of the 
CP in binary lattices was investigated. The expected dy- 
namical power-law scaling relations [see Eqs. (I2~4"|) ] were 
verified and used to obtain the resulting critical expo- 
nents for several sets of parameters (ea, Eb) on the phase- 
separation lines for all the lattices. The evaluation of 
the exponents was performed following Ref. [1] through 
extensive numerical simulations performing averages of 
3 x 10 6 — 10 7 runs to a maximum of t = 3000 time steps. 
Our results obtained for the different lattices indicate 
that, within error bars, the exponents in all cases coin- 
cide with those established for 2d processes in the DP 
universality class (rj = 0.2295(10), 6 = 0.4505(10),2/z = 
1.1325(10)) (25[. Furthermore, the static scaling expo- 
nent ratios determined analogously to the disordered case 
from the QS simulation method are found to coincide 
with those of the DP universality class. 



V. DISCUSSION 

Looking back at the phase diagrams for both the disor- 
dered and the periodic systems, the introduction of dis- 
order in the form of a random placement of A and B 
sites appears to enhance the activity of the system. In 
Fig. [TJ MC data for the disordered system are presented 
and one notices a shallow initial increase in critical Eb 
for a given ea as one moves away from the homogeneous 
critical point followed by an increasingly steep increase at 
values of ea 1^ 0.5 (ea < 0.3). Comparing this behavior 
with the corresponding periodic system [Fig. 3(b)| , the 
critical value for eb in the disordered system is found to 
be much larger. 

Considering the arrangement of, say, A sites as a site 
percolation problem, one notices that for concentrations 
below the percolation threshold (x c ~ 0.59) no infinite 
cluster of such sites can exist. Therefore, no matter how 
small (but non-zero) the corresponding recovery rate ea 
is, it will require a finite value of Eb to render the system 
critical as a finite cluster cannot support an active state 
indefinitely. Conversely, above the percolation threshold 
there exists a finite value of ea below which the system 
will be active irrespective of the value of Eb- There- 
fore, for the case of x = 0.5, no asymptote at any non- 
zero value of ea would be expected. Interestingly, the 
mean field expression Eq. ([7]) does predict an asymptote 
at ea — x albeit for any concentration of sites. 

Turning to the CP in heterogeneous periodic lattices, 
for a range of cases the combination of the standard mean 
field approximation and our alternative analytical ap- 
proximation is useful in practice to pinpoint the location 
of the transition a priori. Indeed, a tight fit for all cases 



(with the exception of the simple chessboard lattice) can 
be attested. In particular, the influence of spatial struc- 
ture on the quality of our approximation becomes evi- 
dent. The less mixed the arrangement of A and B sites 
becomes, the worse the fit of the approximation is found 
to be as indicated by the results in Fig. |3(f)| In order 
to evaluate the practical relevance of our approximation, 
the question of whether it is expected to always yield a 
lower bound has to be addressed. To this end, we define 
an average clustering coefficient specific to a particular 
lattice configuration and site type. For sites of type A for 
instance, define Ca = ^nn,a/^ where rnn.a denotes the 
number of nearest neighbors of site type A (and analo- 
gously for B sites) . For the case of a periodic lattice with 
a 1 : 1 mixture of A and B sites, consider the minimally- 
clustered configuration, that is the standard chessboard 
[lattice (i)], for which Ca = Cb = 0. We know that 
for this case our approximation yields a very tight lower 
bound to the true curve of critical points. Any lattice 
with the same concentration of sites will necessarily have 
a higher clustering coefficient, i.e. a larger fraction of 
contiguous regions of A and B sites. Assuming different 
recovery rates for the two types of site, the disease will 
have the tendency to survive longer in a constellation 
AABB as compared to ABAB due to the adjacency of 
two sites of lower recovery rate (say A) which enhances 
the probability of infection and reinfection in the ^4^4 ar- 
rangement. This activity-enhancing effect is not offset 
by the fact that two less reactive sites (say B) are also 
bordering as their faster (than A sites) recovery is largely 
independent of spatial arrangement. Indeed, a direct di- 
agonalization of the corresponding Liouville operator for 
these two different arrangements of sites readily confirms 
this intuition. One obtains a lower (absolute) value for 
the real part of the first non-trivial eigenvalue in case of 
an arrangement AABB as compared to ABAB indicat- 
ing a slower approach to the absorbing state in a finite 
system. 

From this we conclude that for any periodic arrange- 
ment of A and B sites our alternative analytical approx- 
imation is expected to yield a lower bound to the phase 
separation line. Similarly, the arrangement used in lat- 
tice (iv) is the minimally-clustered (Ca = 1/2, Cb = 0) 
arrangement with a 1 : 2 concentration of sites and is 
found to give a lower bound leading us to expect the 
same behavior for any arrangement of A and B sites in 
this ratio. Therefore, by testing the minimally-clustered 
case for the desired concentration one should in practice 
be able to verify whether or not a lower bound is expected 
by our approximation. 

Considering the critical exponents obtained for both 
the disordered and the periodic systems, they are found 
to be disorder dependent in the former case while they 
remain at their DP values in the latter. Predictions 
from a numerical implementation of the strong-disorder 
renormalization scheme in 2d p redict an exponent value 

•^stroni 



1.0 at an IRFP 



id p redn 



At the same time, 



as conjectured in [19j and supported by numerical evi- 
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dence, the DCP in 2d is likely to be dominated by such 
a fixed point for sufficiently strong disorder similar to 
the Id case. Thus, one would expect disorder-dependent 
varying exponents which approach their values expected 
at an IRFP for strong disorder. Our results support this 
picture: we find continuously varying disorder-dependent 
exponents and a finite-size scaling exponent x which is 
compatible with the predicted value at an IRFP for the 
strongest disorder under consideration (e^ = 0.3). 

Regarding the unchanged exponents in the case of pe- 
riodic systems, these findings confirm theoretical argu- 
ments [27] which make a prediction about the universal 
behavior of the CP in heterogeneous and disordered sys- 
tems. Under coarse-graining the heterogeneity present in 
systems such as the heterogeneous periodic lattices con- 
sidered in this paper will eventually become homogeneous 
after a finite number of iterations of the coarse-graining 
procedure. Thus, one would expect the critical behavior 
of the CP to be governed by the conventional clean fixed 
point of the renormalization group transformations Q . 

VI. CONCLUSION 

In conclusion, we have investigated the contact process 
in both heterogeneous disordered and periodic 2d systems 
(binary lattices). The phase diagram has been obtained 
via extensive Monte Carlo simulation. Furthermore, two 
approximations have been successfully used in order to 
constrain a region of phase space which contains the line 
of critical points. First, the mean field approximation 
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